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ABSTRACT 


Optimal control theory suggests maintaining an orbital altitude band for Low- 
Earth-Orbiting (LEO) satellites using periodic thrusting than forced Keplerian motion, 1.e. 
a trajectory obtained by thrust-drag cancellation. Designing guidance algorithms for orbit 
maintenance is complicated by the nonlinearities associated with orbital motion. An 
algorithm developed previously using thrusters firing significantly off the direction of 
motion successfully maintains an orbital band, but is very inefficient. This thesis develops 
two different control strategies based on the osculating orbital parameters, taking a 
conservative approach to keeping within altitude limitations. Thrust is in the local 
horizontal plane, along the direction of flight. Single and dual burn maneuvers are 
considered for various bandwidths and thruster sizes. The dual burn strategy is somewhat 
close to a Hohmann transfer. The specified orbital band is generally maintained, with 
some cases slightly exceeding the upper limit. Propellant consumption for both maneuvers 
is significantly better than previous methods. This thesis shows that forward firing 
thrusters can be used with osculating orbital parameters to obtain efficiencies within 


forced Keplerian motion values. Accesion For 
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I. INTRODUCTION 


The rising cost of placing a satellite in orbit 1s creating a higher emphasis on 
minimization of non-payload mass to facilitate more payload items, and thus capability, on 
a single vehicle. Reducing the mass of propellant required to maintain a desired orbit has 
the additional benefits for the program in that it also reduces the size and mass of the tanks 
needed to store it. Spacecraft configuration, orbit requirements and thruster control logic 
dictate propellant consumption. The propellant mass required for orbital maintenance 1S 
significantly higher for a Low-Earth Orbiting (LEO) satellite than a geosynchronous one 
because of the significant effects of atmospheric drag on the orbit. 

Consider the problem of maintaining a spacecraft within a prescribed orbital band 
in LEO (see Figure 1). The conceptually simplest thruster control logic is that of a Forced 
Keplerian Trajectory (FKT), where the thrust cancels drag. While this method is easy to 
visualize it is difficult to implement because of technological restrictions on drag 
estimation, thrust vectoring and thrust magnitude adjustments. Although Ross and 
Melton, using optimal control theory, show that an FKT is not an optimal maneuver for 


propellant consumption [Ref. 1], it is a good baseline to measure other methods. 





Rmax 
Figure 1. Orbit Radial Bandwidth 
The Hohmann transfer between the lower and upper radial bands is another 


method of thruster control. Two thruster burn sequences are used in this scheme. The 











change in velocity imparted at the lower altitude in the first sequence places the satellite in 
an elliptic transfer orbit that has an apogee of the maximum radius, while the one at the 
upper altitude circularizes the orbit. Both thrust sequences are in the direction of flight. 
Gottlieb's Eccentricity Intercept Targeting and Guidance (EITAG) [Ref. 2] 
program computes thruster burn duration and timing based on the relationship between 
eccentricity, e, and semi-major axis, a, as shown in Figure 2. It integrates forward from 
the current e-a pair and backward from the desired e-a pair. A two burn routine reboosts 
the satellite to the desired circular orbit. The first burn achieves the intersection point of 
the two curves. The vehicle coasts until it reaches the required position in this 


intermediate orbit that enables the second burn to move the satellite along the final curve. 


Second Burn 





Figure 2. Eccentricity Versus Semi-major Axis 

In order to maintain a radial band Pauls [Ref. 3] and Wilsey [Ref. 4] developed a 
single burn control logic based on orbital radius, specific energy and a thruster cant angle 
(the angle the thrust vector makes with respect to the local horizontal). While the method 
successfully controls the radial band it causes the satellite to consume at least three times 
the propellant of the FKT benchmark. 

Since orbital parameters are derived from the radius and its changes, in this thesis a 
control strategy is developed using the osculating (instantaneous) perigee and apogee radii 
for a single burn with the thruster firing in the local horizontal plane along the direction of 
motion. A dual burn logic, using the single burn strategy as the first burn and the flight 


path angle as the control variable for the second burn, is also explored. 
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Ui. FORMULA DEVELOPMENT 


A. DEVELOPMENT OF THE EQUATIONS OF MOTION 

Assuming the satellite is a non-lifting (blunt) body simplifies analysis in that drag ts 
the only uncontrollable nonconservative force being considered. Specifying planar 
motion, the initial orbit as circular, defining the problem as a two-body problem and 
neglecting all other possible perturbations facilitates ease of formulation. The orbital 


dynamics can be considered in a polar coordinate system as shown in Figure 3. 





Figure 3. Orbital Coordinate System 


The equations of motion are 
F, 
ar=Liz (1) 
F 
at= mo = | (2) 


where a, and a, are the radial and transverse components of the inertial acceleration, LF, 
and =F are the sums of the respective forces, and m is the satellite mass. Drag and thrust 


can be broken into components 





D,;=—-Dcosy (3b) 
T7,/=Tsina (4a) 
7;=Tcosa (4b) 


where y is the flight path angle (the angle between the transverse axis and the velocity 
vector) and a is the thruster elevation angle (the angle between the transverse axis and the 
thrust vector). Expressing the equations of motion in polar coordinates and substituting in 
Equations 3 and 4 yields 


ss oe Dy... a 
r—0 pa t+ 345 (5) 


Ort+2 Or +2 (6) 
where u is the geocentric gravitational constant and @ is the angular position of the 
satellite as depicted in Figure 3. 
B. NONDIMENSIONALIZATION OF THE EQUATIONS OF MOTION 

It is useful to nondimensionalize the equations of motion listed above in order to 
limit the effects of computational errors when dealing with large numbers. This also 
allows better analysis of the effects of the variation of parameters. 

1. Definitions 

Base units in mass, length and time are required to nondimensionalize the problem. 
The base mass, m,, is the initial mass of the satellite. The base length, 7,, is the initial 


semi-major axis of the actual orbit. Since the satellite travels above and below the center 





of the desired radial band the initial satellite position is at the upper limit of the bandwidth. 


The base time, ¢,, is the period of the initial orbit 


3 
” "b 
i nz (7) 


The nondimensionalized mass, radius and time are thus 


— m 
M = hh, (8) 
FS - (9) 
t=- 10 
-. tp ( ) 
a Equation Nondimensionalization 


The satellite's position variables and their derivatives are nondimensionalized using 


the base units 











pa Ga Oe = (11) 

e=86 (13) 

g= @ = Oe BHT (14) 

o= = (2) = 308) = tat) = 2S = 78" 9 





where primes represent a differentiation with respect to nondimensionalized time. 


Substituting these equations and Equations 3 and 4 into the equations of motion yields 





2 | | 
rept _ (15!) - 4» Dsiny , Tsina 
ar 0 hire Cnt mi +o (16) 
WAM so lplross — Posy , Toosa 
20 rpr +20 ae eT (17) 


Rearranging terms 








2 ) | 
wu _% Dsiny ts Tsino 
2 rhbMb mM rhmM, m@ (18) 





5” = 78% t Dcosy t T cosa 








> 0AM Sor rhM, WF (19) 
Substituting in Equation 7 for the base time the equations become 
=f | ie ta) . ts Desiny t, Tsine 
PaO TONS) ~ remem rome 
=~// ‘of t?7 Decos fe 
O oe, ieee b uf » Tcosa (21) 


r rom, mr roMb mr 


3. Nondimensionalization of Nonconservative Forces 
The nonconservative forces (thrust and drag) are to be nondimensionalized by 
multiplying and dividing by the base units required to remove the dimensions. 
a. Nondimensionalized Thrust 


Nondimensionalized thrust is defined by 





2 
lp 


T=s35T (22) 





b. Nondimensionalized Drag 


Nondimensionalized drag is given as 





== f 
D=7%;D (23) 


4, Nondimensionalization of Other Parameters 
The drag the satellite encounters depends solely on the atmosphere the orbit passes 
through. Drag is defined by 


D= +pACpv* (24) 


where p is the atmospheric density, A is the effective surface area of the vehicle, C,, is the 
drag coefficient associated with the spacecraft, and v is the instantaneous velocity. An 


exponentially decaying density is a common atmospheric model and is given by 


1 
p= psc r (25) 


where p, is the density at ,,, the reference altitude, and //is the atmospheric scale height. 


Nondimensionalized scale height is defined as 

= 1% 
p=5 (26) 
A base density is defined as 


Db = Po (27) 











The nondimensionalized density is 


The ballistic coefficient is defined as 


thus the nondimensionalized ballistic coefficient 1s 





Dp. B _ ms 
B= Pols ~~ ACaPpors 


The spacecraft's velocity is 


v2 =r* +7? 7 


and is nondimensionalized using Equations 9, 11 and 14 to 


The nondimensionalized thrust-to-drag ratio 1s 


p= 


where D, is the nondimensionalized drag for the FKT orbit 
“By 





arin 


BH 1 
_ P _ Pe ralt—rre) 2 P 
P=, = po® ae w =e 


(28) 


(29) 


(30) 


(31) 


(32) 


(33) 


_ The change in mass is given 





Lae ae (34) 
Nondimensionalizing yields 
aol ny, fo TP Timers ty Try 
ee Ispg Tp? ™b ~ — Isp8to ae 
5. Nondimensionalized Equations 


Substituting Equations 22 through 32 into the equations of motion produces 


2 ee 
_ =f). e "b / H gin Tei 
NW _6 me _ 1 oe 


r = = 36 
r 2mB (36) 
HW 6’ - Ae) v* cos y T 
a) = -—-,—- - oe tt aces’ 37 
- 2mrB mr 0 
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Ill. COMPUTER MODEL DEVELOPMENT 


A. ORIGINAL DEVELOPMENT 

The original program was developed by Pauls [Ref. 3], and modified by Wilsey 
[Ref. 4] and Gardner [Ref. 5]. It was FORTRAN based, and simulated the motion of a 
spacecraft in orbit by utilizing a fourth-order Runge-Kutta numerical integration routine. 
Wilsey [Ref. 4] nondimensionalized the equations and used four state variables to reduce 
the two second order equations of motion to four first order equations, Gardner [Ref. 5] 
redefined the variables and added the mass to the state variables and the mass equation as 
a first order differential integrated by the Runge-Kutta routine. 
B. STATE VARIABLE DEFINITIONS 

Maintaining a four element state vector based on the position variables and leaving 


the mass as a separate issue, the state vector is defined as 


A eg (38) 
X72 =r=X}] (39) 
x3=0 (40) 
x4 =0=*3 | (41) 


Their nondimensionalized counterparts are 
Xx, =r (42) 
X2=P7 =x; (43) 


x3=6 (44) 








6’ — x" (45) 


C. NONDIMENSIONALIZED STATE VARIABLE EQUATIONS OF 
MOTION 
Rewriting the equations of motion (Equations 38 and 39) by substituting in the 


nondimensionalized state variables, Equations 42 through 45, yields 


2 (2) | ~ 
=f _= a2 _[2n)°_e < ?*¥siny | Tsing 
X54 = X1X4 x mB mn (46) 
a) : 
—/ 2X4x2 © | : | cos T cosa (47) 
4 Xj 2mx,B mx 


The flight path angle, y, is dependent on the radial and tangential components of the 


velocity vector and can be expressed in terms of the state variables 





siny = > (48) 
Xix 
cosy => (49) 
Thus the equations appear in the computer program as 
= 
xj =X? (50) 
—~2 4n? e (-z) wx. . Tsna 

Y= 1X, oo oOo 51 

eed 2mB mi a 
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x5 =X4 (52) 


2X 4X ee) vx T 
4X2 4 cos &@ 
nr ee (53) 
D. PROGRAM ORGANIZATION 

The program, written in MATLAB, is broken into several functional subroutines to 
provide a fluid and logical flow of data through required calculations. Initialization of 
constants and variables is followed by the computation of drag at the Keplerian orbit 
radius and then the main program. Nondimensionalized velocity is computed in LEOVEL 
using the radial and tangential velocity components. The exponential term of drag is 
evaluated in LEODRAG. LEOLOCAL determines the individual pieces of the equations 
of motion for easy assimilation in LEOXDOT, where the derivatives of the 
nondimensionalized state variables are calculated. The Runge-Kutta is completed in 
LEORK4A, LEORK4B, LEORK4C, and LEORK4D, with the program recycling through 
from LEOVEL. After the RK4 routines have produced the new nondimensionalized state 
variables the propellant consumption for that step is computed, then the osculating orbital 
parameters are determined by LEOPARM1 and LEOPARM2. The thruster firing logic 
routine in effect is called (either LEOTFLPH alone or in combination with LEOTFLAP 
for the dual burn scenario) to determine whether or not the thruster is active during the 
upcoming step. A flowchart for the program is shown in Figure 4. 
E. PROGRAM VALIDATION 

All aspects of the code were verified since the program was written in a different 
language from the previous versions. The process Wilsey [Ref. 4] utilized was duplicated. 


1. Initial Validation with No External Forces 


With no atmospheric affects to hinder orbital motion the specific energy and. 


angular momentum for both elliptic and circular orbits must remain constant. Results for 
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these values are shown in Figures 5 and 6 for the elliptic orbit, and Figures 7 and 8 for the 
circular orbit. Small fluctuations are impossible to see on the plots, so closer analysis 1s 
called for. Comparison of the output of the program to the initial values for these 


arguments, depicted in Figures 9 and 10 for the elliptic case and Figures 11 and 12 for the 


Initialize Determine] | Determine 

Variables and [7 Keplerian Velocity 
Determine 
Drag 
Determine 
Differentials 


RK4 
Done? 


Yes 
Propellant 
Consumption 


Figure 4. Program Flowchart 


Constants Drag 








circular one, shows minor fluctuations on the order of 10° with respect to the base for the 
elliptic orbit and zero for the circular orbit. The fluctuations follow no discernible pattern 
and are of a relative magnitude that is known to be beneath MATLAB's computational 
limits and are viewed by the processor as near zero. This suggests they are computational 
noise generated by the numerous calculations of infinitesimal numbers, and the main body 


of the program is validated. 
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Figure 5. Specific Energy of a Drag Free Elliptic Orbit 
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Figure 6. Angular Momentum of a Drag Free Elliptic Orbit 
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Figure 7. Specific Energy of a Drag Free Circular Orbit 
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Figure 8. Angular Momentum of a Drag Free Circular Orbit 
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Figure 10. Change in Angular Momentum of a Drag Free Elliptic Orbit 
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Figure 11. Change in Specific Energy of a Drag Free Circular Orbit 
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Figure 12. Change in Angular Momentum of a Drag Free Circular Orbit 
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2 Validation of Drag 

The program was run using a constant atmospheric density model. The change in 
semi-major axis and velocity were compared with analytic values obtained by using 
equations from Wertz [Ref. 6] that approximate the changes for a satellite in a circular 


orbit experiencing drag 
AC 
Aa = -2n=“pa’ (54) 
AC 
Av = T—pav (55) 


The equations were modified to facilitate ease of computation by making use of the 
definitions of the ballistic coefficient and its nondimensionalized counterpart, Equations 29 
and 30, and noting that the mass and density are constant and that the initial semi-major 


axis is the base radius, to produce 


ACgq 2 m B 
Aa = 21H Poa ae 


a 
BACaBpya B 





(56) 


mp B Vv 
= = T= 57 
BACg Bppa B , ( ) 








_ _ACg 
Av = Tn, PbaVv 


A ten orbit test produced the computational values, the analytic values were calculated 
based on the computational semi-major axis and velocity at the start of each orbit. Percent 
differences were computed for both the change and the element. The results of the test, 
shown in Table 1, show that the percent difference in all cases is less than six-tenths of one 


percent, and the drag routine is validated. 
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Table 1. Percent Difference in Change and Values for SMA and Velocity 
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IV. DEVELOPMENT OF THE ORBIT MAINTENANCE PROCEDURE 


Ross and Melton [Ref. 1] demonstrated the non-optimality of a Forced Keplerian 
Trajectory with respect to propellant consumption. Pauls (Ref. 3] and Wilsey [Ref. 4] 
showed the adequacy, but non-optimality of a single, off axis burn controlled by radius and 
specific energy. A different method of control must be considered in order to approach 
optimality while maintaining the required orbital band. A review of orbital mechanics and 
Wilsey's [Ref. 4] plots for radius and eccentricity suggested a potential solution. 

A. THRUSTER CANT ANGLE 

A firing thruster imparts energy to the satellite at the expense of propellant, thus 

maximizing the energy gained while minimizing the burn time will reduce propellant 


consumption. The energy equation for an orbital body is 


(58) 


The second term of the equation dominates, and orbital energy is negative. As the radius 
decreases the second term grows, and the energy increases negatively. Increasing the 
radius decreases the magnitude of the energy. Zeleny [Ref. 7] notes that Equation 58 can 


be expressed as 


> Sok 
E=>VeV—--Fz (59) 
The change in specific energy with respect to time ts 
de > dv, bar_ L dr 
= OV OG tag = vacos (=o) 27 (60) 


2] 








where & is the angle between the acceleration vector and the tangential axis. Thrust is 
larger than drag because a thrust equal to drag would require a FKT to maintain orbit. 
Figure 13 shows the results at two different thruster cant angles. If the cant angle is small, 


the velocity will not move radically above the tangential axis, and the second term of 





Figure 13. Resulting Velocities at Two Thruster Cant Angles 


Equation 60 will be very small in comparison with the first term. If the cant angle is high 
the initial response will be reduced, since both terms are small. A sufficiently long burn 
would place the velocity and the thrust almost in the same direction, but the objective is to 
reduce burn times. The maximum change of orbital energy is obtained when the velocity 
and acceleration are nearly parallel. For decay of a circular orbit due to drag the flight 
path angle remains small, between 0° and -1°, provided the initial altitude, bandwidth 
(separation between maximum and miniumum radii) and vehicle configuration are 
sufficient to preclude extreme affects. If the path becomes elliptic, or the thruster is small 
and the bandwidth larger than the thruster can achieve in one orbit, the flight path angle 
may vary between 1° and -1°. Setting the thruster cant angle, a, equal to zero enables the 
satellite to sufficiently handle both positive and negative flight path angles with equal 


ability, more so than setting it 1° to either side. To achieve the same change in energy a 
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thruster firing at a significantly larger angle would have to fire longer, consuming more 
propellant. The thruster cant angle is therefore set at zero. 
B. ECCENTRICITY CONTROL 

The imposition of a maximum and minimum radius limits the eccentricity the orbit 
may obtain without having a portion of the orbit exceed either boundary. If the 
eccentricity surpasses the limit, the thruster control logic, whatever it may be, will be 
forced to fire the thrusters more often than if the eccentricity remains within orbital 


bandwidth requirements. The limiting eccentricity is given by 


e _ Rmax-Rmin 
max “~ Rmax+Rmin 


(61) 


For a thruster control logic to be remotely close to optimal it must avoid allowing the 
osculating eccentricity exceeding the maximum value. 
C. SINGLE BURN THRUST CONTROL LOGIC 

Since the orbital parameters of interest are based on the radius, orbital angle and 
their respective rates of change, and since the desired outcome is the maintenance of a 
prespecified orbital band by limiting the eccentricity, the control logic utilizes these values 
to obtain the desired results. 

1. Thruster Firing Criteria 

Acknowledging that the satellite would descend beneath the minimum radius if the 
thruster was not fired prior to reaching it dictates a burn start before the lower radial limit 
is reached. Atmospheric density profile, orbital bandwidth, thruster capacity and cant 
angle, and satellite configuration all play a role in how far the vehicle drops from thruster 
initiation until the radius ceases to decrease. Starting the burn when the osculating perigee 
radius drops below the minimum radius will at worst delay the time until the satellite 
passes the minimum value or cause the satellite to begin to climb well before the lower 
radial limit is reached. Any result in between the two is reasonably acceptable and the 


thruster firing criteria is established. 
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Z; Thruster Turn Off Criteria 

Gottlieb [Ref 2] shows the relationship between semi-major axis, @, and 
eccentricity, e, during thrusted flight. Ifthe thruster is allowed to fire until either portion 
of the a-e pair is beyond its limiting value the satellite will in all likelihood either exceed 
the bandwidth restrictions or require excessive burns to maintain it. Turning the thruster 
off when the osculating apogee radius exceeds the maximum radius implies that the orbital 
radius will never exceed the maximum value. This is mainly because the effect drag will 
have on reducing the apogee radius, and partly due to the change in apogee for any given 
instant being small enough so that the osculating apogee does not significantly exceed the 
limit. The error in overshooting the upper limit is small considering the second factor, and 
the thruster turn off criteria is established. Figure 14 is the logic flowchart for the single 
burn maneuver. 

3. Resulting Orbit 

Since the thrusts occur based on the osculating perigee falling below the minmum 
limit the satellite may not be at a perigee position when the thruster 1s activated. 
Additionally, the velocity change will not necessarily correspond to that of a Hohmann 
transfer. Changes in the perigee radius and apogee radius may occur simultaneously. 
Raising the instantaneous radius by increasing the apogee radius is the primary intent of 
the burn, but changes in the perigee radius are also important. Raising the perigee radius 
along with the apogee radius places the satellite in a more desirable orbit because it 
increases the time it takes for the perigee radius to decay below the minimum limtt, 
increasing the period between burns. 

The change in perigee radius can be determined from analysis of the semi-major 
axis, eccentricity, angular momentum and specific energy. Zeleny (Ref. 7] states the 


change in eccentricity and angular momentum from a change in velocity are 


de = ade + 2ehdh) (62) 


24 








— 
dh= r xav (63) 
where A is the angular momentum. The angular momentum is 


7 > > 
h = v 


r x (64) 
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Figure 14. Single Burn Maneuver Thruster Control Flow 
and thus 
h=rvsin(2-y) = rv COSY (65) 


The change in angular momentum can be expressed as 
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dh = racosé (66) 
Substituting Equations 58, 60, 65 and 66 into Equation 62 yields 
2 
Y 
de = -alrv'a cosy cos (y — &) + “cosy 


+rvacos & — 2ua cos €] (67) 


Since all the angles in this equation are zero or approximately zero, and the terms outside 


the brackets are positive, the sign of the equation can be determined from 
. . ep) uv? 
sign(de) =sign| 2a(rv* — W) + (68) 


Zeleny [Ref. 7] states that the change in semi-major axis is 


da = = de (69) 
The perigee radius is defined as 
rp = a(1—e) (70) 
so its change 1s 
dr, = (1 — e)da — ade (71) 


From Equation 68 the change in eccentricity varies from positive to negative. As the 


velocity increases with a burn, the change in eccentricity is more likely to become positive, 
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and attempt to reduce the change in perigee radius. Eccentricity is generally very small, 
and changes to it even smaller for any given instant in time. The change in semi-major 
axis is not small, in comparison. The first term in Equation 71 is on the order of the 
change in semi-major axis while the second term is significantly smaller. The overall result 
is a positive change in the perigee radius, placing the satellite in a much more desirable, 
but still eccentric, orbit. 
D. DUAL BURN CONTROL LOGIC 

The single burn strategy described above leaves the satellite in an elliptic orbit. If 
the perigee radius has not been sufficiently raised, drag will soon lower it beneath the 
minimum limit and the thruster will fire again. Reducing the eccentricity by conducting a 
second burn as the satellite approaches apogee could result in an increase in the period of 
radius boosting burns, potentially saving propellant. The objective of the second burn is to 
decrease eccentricity, but not necessarily make the orbit circular. 

l. Thruster Firing Criteria 

As the satellite approaches apogee the thruster must fire to increase both the 
semi-major axis and the perigee radius. The flight path angle, y, 1s positive during the 
perigee to apogee transition, and decreases as the vehicle nears apogee. Waiting for the 
flight path angle to become negative before initiating thrust could raise the apogee radius, 
placing the satellite on a path that exceeds the maximum radial limit. Firing the thruster 
when it drops below a minimum positive value seems to be the prudent course of action. 
An arbitrary value was selected for the firing criteria, and adjusted until it minimized 
eccentricity for the smallest thruster-bandwidth combination. Using the resulting firing 
angle with larger thrusters and bands yielded higher than desired eccentricities. More 
powerful thrusters change the path quicker, and should require less time to reduce the 
eccentricity. Larger bandwidths increase the maximum allowable eccentricity, producing 
boost orbits with a correspondingly lower velocity at apogee. Longer burns are required 
at the smaller velocity to raise it to a value close to that of a circular orbit at the apogee 


radius. Modifying the base flight path firing angle by a ratio of the bandwidth-to-thrust 
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ratio proved insufficient. Using the square of the bandwidth in the numerator of the ratio 
produced more favorable, but not perfect, results. 

z Thruster Turn Off Criteria 

If the thruster fires too long the apogee radius will rise, and the path potentially 
exceed the upper radial limits. If the burn is too short the perigee radius will not be 
increased sufficiently to affect its decay-to-minimum time. Controlling the end of the burn 
by comparing the perigee radius to the radius is insufficient because it does not consider 
the effects on the apogee radius, and the burn might be long. Examination of the flight 
path path angle proves more interesting. 


The tangent of the flight path angle ts given by 


= eee os 
¢ = tany= 3 (72) 
Its rate of change is 
dr dr rd0 (i 2 i) 
eS ee pe r, 
ds r@ r*6 -@? ~~ FO —? 9 W?) 


Near apogee, with no thrust, the values of the variables inside the expression are 


r>>0 (74a) 
r—> 0" (74b) 
r< 0 (74c) 
Q> 0 (74d) 
6< 0 (74e) 
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and the change in the tangent of the flight path angle is negative. The thruster, as seen 
from the equations of motion, changes the rates in Equations 74c and 74e to positive 
values. Since radial velocity is approximately zero, the radial acceleration changes govern 
Equation 73, eventually causing d¢ to become positive. 

Continuous thrusting at the chosen cant angle increases both r and v. At apogee 
thrusting has more affect on velocity than radius. The result is that the orbit tends to 
circularize, then, as shown previously, the semi-major axis increases, and the orbit 
becomes elliptic again. This ties in with the change the flight path angle undergoes, and 
thus when it begins to increase, the eccentricity has reached a minimum. The thruster cut 
off criteria is determined. Figure 15 depicts the flowpath for the second burn in the dual 
burn thruster control logic. 

3. Resulting Orbit 

The resulting orbit will be nearly circular and at or near the maximum radial limit. 
Since the perigee radius is almost equal to the apogee radius, the decay pattern should 


allow for more time between boosting burns. 
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Figure 15. Dual Burn Strategy Second Burn Logic 
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Vv. ANALYSIS AND RESULTS 


The satellite is initially positioned at the upper radius limit, as though a launch 
vehicle had placed it there. FKT analysis is done at the lower limit, midband, and the 
initial radius. Configuration changes are limited to the thruster sizing and the limited 
effects on nondimensionalization caused by the different starting altitudes. The following 
parameters were used in all cases: 


Midband Radius: R, + 260 km 


Cant Angle: 0° 

Ballistic Coefficient: 150 kg/m’ 

Scale Height: 46.9 km 

Specific Impulse: 300 s 

Initial Mass: 20,000 kg 

Base Density Altitude:250 km 

Base Density: 7.248 x 10" kg/m” 
Density Factor: 12.47 


Bandwidths were selected as 2, 10, 25, and 75 km (0.3 x 10°, 1.5 x 10°, 3.8 x 10°, and 
11.2 x 10° distance units, respectively). Thruster sizes were limited to 40, 80, 160, 320, 
640, and 1280 N. Table 2 shows the nondimensionalized counterparts to thruster size. 
The program is run for 100 orbits for the first three bandwidths, and 200 orbits for the last 
bandwidth. The equations of motion are updated 1000 times and the output sampled ten 
times per orbit. 
A. SINGLE BURN STRATEGY 

1. Band Maintenance 

Pauls [Ref. 3] and Wilsey [Ref. 4] demonstrated that their single burn control logic 
did not permit forward-firing thrusters to maintain an orbital band. The ability of the 
selected Thruster Firing Control Logic (TFCL) to keep the satellite in a bandwidth at all 
must first, therefore, be established. Initial runs demonstrated that the logic successfully 


maintained the radius between two boundaries. Having shown the capacity to constrain 


3] 





the position, the TFCL had to be evaluated on its ability to hold specified limits. 
Variations in bandwidth and thrust magnitude were applied. Figures 16 through 39 show 
the radius versus orbit. The TFCL kept the vehicle at or above the minimum radius 
regardless of the band and thruster size selected. Using thrusters smaller than 40 N could 
cause the path to drop under the limit. Staying beneath the upper radial limit depended on 
the thrust and the bandwidth, as well as the time step size. High thrust at small 
bandwidths caused the radius to exceed the maximum boundary. Table 3 summarizes the 


single burn band maintenance results. 





Thrust [N] 


es 
se _ 20] _10] _o] _stol_zto 


Bandwidth Nondimensionalized Thrust 
[km] 


a eal iad] _340| 6985] __ 139.69] 27939 
io 874[1748|__3497] 6993] 13987] 270.74 
asl ate[ 17a] 35.05) 7009) 140.18) 28037 
asf gal ates] 3531] 7062) _tai.2q] 242.48 


Table 2. Nondimentionalized Thrust by Bandwidth 
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Table 3. Single Burn Strategy Band Maintenance Summary 
Z. Orbital Path 
The trajectory caused by the single burn, as expected, is generally elliptic. The 
orbital radius follows no immediately discernible pattern, and neither does the burn 
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pattern. In some instances TFCL causes the thruster to fire as the vehicle approaches 
apogee, inadvertently reducing the eccentricity of the orbit, as demonstrated in Figures 40 
and 41. This is because the routine is only looking at the osculating perigee radius for 
firing criteria, not the vehicle's position relative to apogee or perigee. Since the orbit is 
not permitted to exceed the maximum eccentricity imposed by the radial limits, burns are 
required less frequently than by the Pauls-Wilsey method. 

The method utilizes the full bandwidth during the course of the run, but does not 
use all available distance units on every opportunity. In the smaller bands, particularly the 
2 km (0.3 x 10° distance units), the logic uses as little as half the available bandwidth. In 
cases with a larger minimum-maximum separation, approximately 90% of the band is used 
when burns occur. 

B. DUAL BURN STRATEGY 

L Band Maintenance 

Again the logic had to be tested to ensure it was capable of maintaining a band. 
Runs showed the dual burn maneuver kept the vehicle within radial limits. The new 
strategy was put through the same bandwidth-thruster variation sequence the single burn 
control method was. Plots of the resulting radius versus orbit are shown in Figures 42 
through 65, Since the single burn strategy was used as the initial burn, the radius was not 
expected to, and did not, go below the minimum. The timing of the second burn did not 
always meet the needs of the particular osculating orbital elements, and periodically led to 
the upper radial limit being exceeded. Review of the burn initiation time, the flight path 
angle and the satellite's position with respect to apogee revealed thrusting began earlier 
than required. Thrust initiation for the second burn is based on a scaled flight path angle. 
The scaling ratio was too large, indicating the method favors the bandwidth over the 
thrust capacity more than it should. The use of a ratio was an arbitrary choice, other 
functions could well be the solution to the problem. Table 4 summarizes the dual burn 


maneuver results. 
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Figure 16. Single Burn Strategy Orbital Radius (2 km Bandwidth, 40 N Thrust) 
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Figure 17. Single Burn Strategy Orbital Radius (2 km Bandwidth, 80 N Thrust) 
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Figure 18. Single Burn Strategy Orbital Radius (2 km Bandwidth, 160 N Thrust) 
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Figure 19 Single Burn Strategy Orbital Radius (2 km Bandwidth, 320 N Thrust) 
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Figure 21. Single Burn Strategy Orbital Radius (2 km Bandwidth, 1280 N Thrust) 
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Figure 22. Single Burn Strategy Orbital Radius (10 km Bandwidth, 40 N Thrust) 
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Figure 23. Single Burn Strategy Orbital Radius (10 km Bandwidth, 80 N Thrust) 
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Figure 24 Single Burn Strategy Orbital Radius (10 km Bandwidth, 160 N Thrust) 
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Figure 25. Single Burn Strategy Orbital Radius (10 km Bandwidth, 320 N Thrust) 
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Figure 26. Singie Burn Strategy Orbital Radius (10 km Bandwidth, 640 N Thrust) 
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Figure 27. Single Burn Strategy Orbital Radius (10 km Bandwidth, 1280 N Thrust) 
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Figure 28. Single Burn Strategy Orbital Radius (25 km Bandwidth, 40 N Thrust) 


6660 


ie cee ee ee hon a See SS 


: 6655 
6650} 
| 


| 6645 | | | | | 
| Radius Hil TTT Wilh | 
| HH HIN ARETE } 

| [km] 6640 : AM AAA 
| HT AN iy 
| | | 


6635 1 TATA THRU IMM 


6630 





-——_— oon oO eo ee eee ee eee awe oe ow ee wee eee 


| 6625-— 
0 20 40 60 80 100 
Orbits 


Figure 29. Single Burn Strategy Orbital Radius (25 km Bandwidth, 80 N Thrust) 
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Figure 31. Single Burn Strategy Orbital Radius (25 km Bandwidth, 320 N Thrust) 
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Figure 32. Single Burn Strategy Orbital Radius (25 km Bandwidth, 640 N Thrust) 
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Figure 33. Single Burn Strategy Orbital Radius (25 km Bandwidth, 1280 N Thrust) 


42 








ie! keene Ne 


Figure 35. Single Burn Strategy Orbital Radius (75 km Bandwidth, 80 N Thrust) 
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Figure 36. Single Burn Strategy Orbital Radius (75 km Bandwidth, 160 N Thrust) 
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Figure 37 Single Burn Strategy Orbital Radius (75 km Bandwidth, 320 N Thrust) 
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Figure 39 Single Burn Strategy Orbital Radius (75 km Bandwidth, 1280 N Thrust) 


45 















Apogee/ 
Radius/ 
Perigee 


[km] 









Figure 40. Apogee. Radius, and Perigee (10 km Bandwidth, 80 N Thrust) 
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Figure 41. Apogee, Radius, and Perigee (25 km Bandwidth, 40 N Thrust) 
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Figure 43. Dual Burn Strategy Orbital Radius (2 km Bandwidth, 80 N Thrust) 
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Figure 44. Dual Burn Strategy Orbital Radius (2 km Bandwidth, 160 N Thrust) 





40 60 80 100 





Orbits 





Figure 45. Dual Burn Strategy Orbital Radius (2 km Bandwidth, 320 N Thrust) 
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Figure 47. Dual Burn Strategy Orbital Radius (2 km Bandwidth, 1280 N Thrust) 
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Figure 49. Dual Burn Strategy Orbital Radius (10 km Bandwidth, 80 N Thrust) 
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Figure 50. Dual Burn Strategy Orbital Radius (10 km Bandwidth, 160 N Thrust) 
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Figure 51. Dual Burn Strategy Orbital Radius (10 km Bandwidth, 320 N Thrust) 
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Figure 53. Dual Burn Strategy Orbital Radius (10 km Bandwidth, 1280 N Thnust) 
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Figure 55. Dual Burn Strategy Orbital Radius (25 km Bandwidth, 80 N Thrust) 
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Figure 56. Dual Burn Strategy Orbital Radius (25 km Bandwidth, 160 N Thrust) 





Figure 57 Dual Burn Strategy Orbital Radius (25 km Bandwidth, 320 N Thrust) 
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Figure 58 Dual Burn Strategy Orbital Radius (25 km Bandwidth, 640 N Thrust) 
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Figure 59. Dual Burn Strategy Orbital Radius (25 km Bandwidth, 1280 N Thrust) 
55 








Figure 61. Dual Burn Strategy Orbital Radius (75 km Bandwidth, 80 N Thrust) 
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Figure 62. Dual Burn Strategy Orbital Radius (75 km Bandwidth, 160 N Thrust) 
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Figure 63. Dual Burn Strategy Orbital Radius (75 km Bandwidth, 320 N Thrust) 
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Figure 65. Dual Burn Strategy Orbital Radius (75 km Bandwidth, 1280 N Thrust) 
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Z. Orbital Path 

The orbital trajectory for the dual burn maneuver was significantly smoother than 
that of the single burn. Raising the perigee radius with the additional thrusting sufficiently 
reduced the eccentricity 5 enable a near-circular decay pattern to evolve in most cases. In 
those where the burn started early the eccentricity was decreased enough to prevent the 


wide range of apogee and perigee experienced in the single burn maneuver. 
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Table 4. Dual Burn Strategy Band Maintenance Summary 


C. PROPELLANT CONSUMPTION 

Propellant consumption is calculated for each variation of the single and dual burn 
maneuvers, as well as for an FKT at the maximum radius, midband radius, and minimum 
radius. An additional consideration involves the strategy of using a FKT to keep the 
satellite at the upper radial limit until the time when removing any thrusting would cause 
the orbit to decay to the minimum radius by the end of the run period. TFCL consumption 
is compared against these four base usage patterns. 

Consumption for both the single and dual burn strategies is approximately the 
same. The exact value at any given instant may be higher for one or the other, but the 
average values overall are generally equal. Plateaus between burns are longer in the dual 


burn case, as expected. Burns occur more frequently in the single burn maneuver because 


59 





the perigee radius may not be sufficiently raised by any given burn to preclude drag 
quickly decaying it below the minimum limit. 

Comparison of the two methods with the four baselines described above provides 
insight into their efficiency. In all cases considered the propellant required by either TFCL 
strategy is slightly less than or roughly equal to the midband FKT counterpart. As the 
bandwidth gets very large the TFCL values appear to drop away from the midband ones. 
The low altitude FKT consumes more, significantly more in cases of larger bandwidths, 
than the midband FKT. The upper radius FKT and its modified version use less that the 
midband FKT. Both TFCL models seem to be bounded above by the midband FKT and 
below by the high altitude FKT, as indicated in Figures 66 through 89. This suggests that 
although TFCL is better than attempting an FKT at the middle of the radial bandwidth, it 
is less efficient than either an FKT at the maximum radius or the modified maneuver 


described previously. 
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Figure 67. Propellant Consumption (2 km Bandwidth, 80 N Thrust) 
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Figure 68. Propellant Consumption (2 km Bandwidth, 160 N Thrust) 
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Figure 69. Propellant Consumption (2 km Bandwidth, 320 N Thrust) 
62 | 





600 


Dash-Dot - Single Burn 
Dotted - Dual Bum 


500 
| Solid - Midband FKT 
| Dashed - Max/Min FKT 
4oo+ (Max FKT also with decay to mm) 
| Propellant 
Consumed 
300 
[kg] 
200 
100 





0 20 40 60 80 100 
| Orbits 


C:ash-Dot - Single Burn 
Dotted - Dual Burn 
Solid - Midband FKT 
Dashed -Max/Min FKT 
(Max FKT also with decay to min) 












500 









400 


Propellant 
Consumed 


[kg} 


Figure 71. Propellant Consumption (2 km Bandwidth, 1280 N Thrust) 
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Figure 73. Propellant Consumption (10 km Bandwidth, 80 N Thrust) 
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Figure 75. Propellant Consumption (10 km Bandwidth, 320 N Thrust) 
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Figure 76. Propellant Consumption (10 km Bandwidth, 640 N Thrust) 
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Figure 77. Propellant Consumption (10 km Bandwidth, 1280 N Thrust) 
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Figure 79. Propellant Consumption (25 km Bandwidth, 80 N Thrust) 
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Figure 81. 
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Figure 83. Propellant Consumption (25 km Bandwidth, 1280 N Thrust) 
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Figure 85. Propellant Consumption (75 km Bandwidth, 80 N Thrust) 
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Figure 86. Propellant Consumption (75 km Bandwidth, 160 N Thrust) 
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Figure 87. Propellant Consumption (75 km Bandwidth, 320 N Thrust) 
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Figure 88. Propellant Consumption (75 km Bandwidth, 640 N Thrust) 
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Figure 89. Propellant Consumption (75 km Bandwidth, 1280 N Thrust) 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


The purpose of this thesis was to develop a new thruster firing control logic for 
either a single burn or dual burn maneuver that successfully maintained an orbital band 
while achieving propellant efficiency approaching or exceeding that of forced Keplerian 
motion. Nondimensionalization of the equations permitted more accurate computation 
and easier visualization in variation of parameters. Results indicate that both strategies are 
roughly as efficient as a midband FKT, but less efficient than an upper limit FKT or an 
upper limit FKT with decay to minimum radius at end of life. Reducing the eccentricity 
with a second burn increases the period between boosting burns, but does not reduce 
propellant consumption. 

The question of optimality of thrust cant angle is opened again, based on the new 
thruster firing control logic. Earlier efforts using the instantaneous radius and osculating 
specific energy indicated that a high angle was required to keep a vehicle within a radial 
band, and that low angles were incapable of constraining the radius. Thrusting along the 
transverse axis using the osculating perigee and apogee radii as initiator and terminator, 
respectively, is capable of orbital band maintenance and significantly more efficient than 
previous methods. 

In a separate but related strategy an additional burn, based on the flight path angle 
as the vehicle approaches apogee, reduces the osculating orbital eccentricity, allowing a 
near-circular decay pattern to develop. The additional thrusting raises the perigee radius, 
enabling a longer time between orbit boosting maneuvers. The single and dual burn 
maneuvers are equally effective at band maintenance and propellant consumption, but the 
second method provides a potentially more desirable orbital pattern. 

It is unclear which FKT model is the appropriate baseline to compare orbital band 
maintenance maneuvers on. At the low altitude of the band the FKT solution consumes 
more propellant because it is fighting a denser atmosphere. Likewise, the high altitude 


case and any derivatives from it have the benefit of having a significantly rarefied medium 
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to travel through, and require less propellant. The midband version is not the most 

efficient, thus its validity as the basis for comparison is also questionable. 
The results of this thesis indicate several areas for future research: 
(1) For apogee rasing burns fire the thruster only when the satellite is at 
perigee. The initial burn in both methods takes place when the osculating perigee 
reaches the minimum radial limit, regardless of the position of the vehicle with 
respect to perigee. Examining at each occurring perigee the radial decay from the 
previous perigee, and firing the thruster if the drop in altitude is larger than the 
remaining distance to the lower altitude limit could reduce propellant consumption. 
(2) Refining the thruster initiation for the second burn. Modification of the 
scaling factor or exploration of other controlling functions could improve the 
performance. 
(3) | Coupling EITAG with a control logic and an orbital propagator. The 
resulting propellant consumption and band maintenance could then be compared 
with FKT, Pauls-Wilsey and this method. 
(4) The effect varying the specific impulse has on band maintenance and 
propellant consumption. 

Finally, additional control strategies not discussed here possibly exist that could provide 


the answer to the optimal control of orbital band maintenance. 
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APPENDIX 


MAIN PROGRAM LISTING 


MAIN PROGRAM 


Variable Definitions 


a 
alphad 
alphar 
alto 

B 

Boar 
beta 

D 

Dbar 
at 

dv 

az 

DK 


ES 

ESO 

Et 

g 
gammad 
gammar 
Gnot 
Gto 

h 

LSD 

L 
lambda 
Ptactor 
m 
mbar 
mf 

mfk 
mit 
mftK 
mnd 

mo 
orbits 
ra 

rho 
rhoalt 
rhofac 
rhonot 
rhostd 
RmaX 
Rmin 


Semi-~Major Axis 

Thrust Angle wrt LH [degrees] 
Thrust Angle wrt LH [radians] 
Midband Altitude [km] 
Ballistic Coefficient [kg/m‘*2] 
Nondimentionalized Ballistic Coefficient 
Atmospheric Scale Height 

Drag [N] 

Nondimentionalized Drag 

Drag Factor 

Change in Velocity [m/s] 
Orbital Bandwidth [km] 

FKM Drag [N] 


= Eccentricity 


i 


Maximum Specific Energy at Maximum Radius [J/kg] 
Minimum Specific Energy at Minimum Radius [J/kg] 
Specific Energy [J/kg] 

Initial Specific Energy [J/kg] 

Total Energy [J] 

Standard Earth Gravity [m/s%*2] 

Flight Path Angle wrt LH [degrees] 

Flight Path Angle wrt LH [radians] 

Base Flight Path Firing Angle [radians] 

Modified Flight Path Firing Angle [radians] 
Angular Momentum [km*2/s} 


= Specific Impulse [s] 


= Length Conversion Factor 


= Mass of Fuel Burned in Time Increment 


Local Computational Variables 

Modified Thrust Firing Scaling Factor 
[1000 m/km] 
Spacecraft Mass [kg] 

Nondimentionalized Mass 

[kg] 


FKM Mass of Fuel Burned in Time Increment [kg] 


= Total Mass of Fuel Burned [kg] 


FKM Total Mass of Fuel Burned [kg] 
Arbitrary Mass Used for Nondimentionalization [kg] 
Initial Spacecraft Mass [kg] 

Number of Orbits Completed 

Apogee Radius [km] 

Calculated Atmospheric Density [kg/m*3] 
Reference Atmospheric Density Altitude [km] 
Density Variation Factor 

Reference Atmospheric Density [kg/m‘*3] 
Standard Earth Density [kg/m‘*3] 

Maximum Radius of Band [km] 

Minimum Radius of Band [km] 
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* Rminb = Nondimentionalized Minimum Radius of Band 


= £6 = Initial Orbital Radius [km] 

* rp = Perigee Radius [km] 

aot = Time [orbits] 

=~ TBW = Base Thrust-Bandwidth Scale Factor 

=x Te = Farth's Surface Rotational Period [s]} 

* CL = Final Step Time [s] 

= St). = Thrust Firing Logic Selector 

= Th = “Thrust. LN 

* Thbar = Nondimentionalized Thrust 

* thetad = Angle from Reference Axis [degrees] 

* thetar = Angle from Reference Axis [radians] 

« ThK = FKM Thrust [N] 

* Thm = Maximum (Blowdown) Thrust [N] 

* Thmbar = Nondimentionalized Maximum (Blowdown) Thrust 
« tinc = Increment of Time (Step Size) [s] 

«OL. = Tollerance Value for Computation 

* Tpo = Initial Orbital Period [s] 

* Tp = Orbital Period [s] 

* tstart = Start Time [s} 

x tstop =-Stop: Time [s] 

= V = Velocity [km/s] 

« Vbar = Nondimentionalized Velocity 

« Vimax = Maximum Velocity at Present Orbit [km/s] 
* Vimin = Minimum Velocity at Present Orbit [km/s] 
+ Vmax = Maximum Velocity at Maximum Radius [km/s] 
= Vmin = Minimum Velocity at Minimum Radius [km/s] 
* Vprev = Previous Velocity [km/s] 

Ss C1) = Orbital Radius [km] 

= «x (2) = Orbital Radial Velocity [km/s] 

* x(3) = Orbital Theta 

* «(4) = Orbital Angular Velocity [1/s] 

* xbar(1) = Nondimentionalized Orbital Radius 

’ xbar(2) = Nondimentionalized Orbital Radial Velocity 
$ xbar(3) = Nondimentionalized Orbital Theta 

* xbar(4) = Nondimentionalized Orbital Angular Velocity 
= xbdot(1) = Time Derivative of xbar(l) 

¢ xbdot(2) = Time Derivative of xbar(2). 

=< xbdot(3) = Time Derivative of xbar(3) 

+ xbdot(4) = Time Derivative of xbar(4) 


clear all 


+» Constants 


mu = 3,98601208133E5; * Earth's Gravitational Parameter 
* [km*3/s%2] 
g = 9.806; * Gravity at Earth's Surface 
= {[m/s%*2} 
Re = 6378.145; * Farth's Radius 
= [km] 
ehmona = pi/l180; =« Change Degrees to Radians 
lfactor = 1000; * Convert km to m 
tol = le-6; 
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,cCOuUNne = 1:4; * State Variable Numbers 


* Independent Variable Initialization 


alto = 260.0; 

az = 25.0; 

a = Retaltotdz/2; 
alphad = 0.0; 

B = 150.0; 

betal = 46,9; 

aft = 130: 

dv = 0.0; 

e = 0.00; 

Isp = 300.07 

index =A 

mo = 20000.0; 

mnd = 20000.0; 
rhoalt = 250.0; 

rhostd = 7,248e-11; 
rhofac = 12.47; 

rhonot = fhofac*rnosta; 
Te = Qe pr 4 GRe "3 ) 7 mu) 0.0) 4 
thetad = 180.0; 

Thm = 320.0; 

Thfact = 1.0; 

eine = 0.001; 

tstart = 0,0; 

tstop = 100; 

c = tstart; 

te = 0; 

Gnot = 20e-6 

TBW = 10; 

lambda = TBW/ (Thm/dz%2) ; 
Gto = lambda*Gnot; 
aan = le-14; 

shoodi = Q; 

doozit = 100; 


2 Initialize Orbital Element Variables 


thetar = thetad*dtor; 

ra = a*(l+e); 

ro = (a*(l-e*2))/(l+te*cos (thetar) ); 
rp = a*(l-e); 

Tpo = 2*pi*sgqrt((a%*3)/mu); 

Vy =. Sart ((2*mu/ ro)}—(mu/a)); 

ao ee ae 

x ti) = ro; 

x (2) = e*(sin(thetar) )*sqrt (mu/ (a* (1-e%2))); 
x (3) = thetar; 

x (4) = (sqre(a*mu*(l-e*2)))7 (e(ly 2) 3 
xbar(1) = x(1)/ao; 


TH 





xbar(2) 
xbar (3) 
xbar (4) 
ESo 

Et 

ES 





HH 


i 


— 





(Tpo*x(2))/ao; 


a a 
(Tpo*x(4)); 
((V*2)/2-mu/ro); 
mo*Eso; 

ESO} 


* Initialize Satellite Variables 


mM 
alphar 
mbar 
Thmbar 
Bbar 
Thbar 
mft 
mftk 


mo}; 
alphad*dtor; 

m/mnd; 
(Thm*Tpo*2) / (mnd*ao*lfactor) ; 
B/ (rhonot*ao*lfactor); * m/km 
04.0% 

O20 

O03 


> Initialize Altitude Band Constraint Variables 


RmaX 
Rmin 
Rminb 
Vmax 
Vmin 
Esmax 
Esmin 
gammar 
norbs (1) 
Vimax (1) 


orad(1) 
oV (1) 
ogni.) 
oh (1) 
oEs (1) 
oe (1) 
oa(l) 
ora ( 
orp { 
oTp ( 
t£im ( 


1) 
1) 
1) 
-) 


iol t Wow tt 


rot+dz/2; 

ro-dz/2; 

Rmin/ao; 

(mu/Rmax)*0.5; 

(muy Rm)? 0.57 

((Vmax*2)/2)- (mu/Rmax) ; 
((Vmin*2)/2)-(mu/Rmin) ; 
atan(xbar(2)/(xbar(1)*xbar(4))); 
0; : 

Vmax; 


/ 


* Set Forced-Keplerian Motion Drag 


rK(1) = (ro-dz/2)/ao; 

vK = Tpo*sgrt ((2*mu/ (rK(1)*ao) )-(mu/(rK(1)*ao)))/a0; 
[rho, Dbar] = leodrag(rK, ao, rhoalt, Re, beta, vK, Bhar; :di}; 
DK = Dbar*mnd*ao*lfactor/(Tpo*%2); * 1000 m/km 
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* Main Program 
while t <= tstop, 
while index <= 4, 


«= Determine Nondimentionalized Velocity 


[Vbar] = leovel (xbar) ; 


* Determine Nondimentionalized Drag 


{rho, Dbar] = leodrag(xbar, ao, rhoalt, Re, beta, Vbar, Bhar, df); 


* Determine Local Variables 

[Ely 2g tS) bad, tS): 2b (6) y.. dehy) = eotocal txbar, “mbar, 
Vbar, Thbar, Dbar, alphar); 

* Determine Positional Differential Values 


[xbdot(1), xbdot(2), xbdot(3), xbdot(4)] = leoxdot(xbar, L); 
~« Perform RK4 


if index == 

[xbt, xbdt, xbart] = leork4a(icount, xbar, xbdot, tinc); 
elseif index == 

[xbdt,. xbart) = Lecrk4b (2 count, -xbt, xbat; xbdot;, tine); 


elseif index == 


[xbdt, xbart] = leork4c(icount, xbt, xbdt, xbdot, tinc); 
else 
[xbart] = leork4d(icount, xbt, xbdt, xbdot, tinc); 
end 
xbar(icount) = xbart(icount); 
if rem(index, 2) == 1 


C= 40 5" cinc 
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end 


index = index+l; 
end 
> Reset index value 


index = 1; 


* Compute Propellant Consumption 


Th = Thbar*mnd*ao*lfactor/ (Tpo*2Z) 3 
mf = "The tine*Tpo0/ (isp*¢G) + 

m = m-mf; 

mft = mft+mf; 


mfK = DK*tinc*Tpo/ (Isp*g); 
mftK = mftK+mfK; 
mbar = m/mnd; 


* Increment Output Counter 


shoodi = shooditl; 


2 Store Propellant Consumption 


if rem(shoodi,doozit) == 0 
mtv(fix(shoodi/doozit)+l) = mft; 
mtKv (fix (shoodi/doozit)+1l) = mftK; 
mrat(fix(shoodi/doozit)+1) = mft/mftK; 


tilmctlis (shododi7dooZit) +1) =trl; 


end 


* Determine Orbit Parameters 
eprev = e; 
grp=gammar; 
Vprev = V; 
[r, V, gammar, h, Es] = leoparml(xbar, Tpo, ao, mu) ; 
[e, a, ra, rp, Tp] = leoparm2(Es, h, mu); 
* Thruster Firing Evaluation 
if tfl == 2 


if waar wep 
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end 
end 
if tlt == 
if Thbar > 0 
if gammar-grp >= 0 
tfl = 1; 
end 
end 
end 
if tfl == 
[Thval, tfl] = leotflph(rp, Rmin, ra, 
elseif rem (tfl,2) == 1 


{(Thval, tfl] = leotflapn(gammar, gIp, 
Gto); 


end 


Thhar = Thval*Thiact; 


1s Oi Oe St a 


cfl = O; 


e = 0; 

end 
+ Output 
if rem(shoodi,doozit) == 0 

orad (fix (shoodi/doozit)+l) = x; 

6V (fix shoodi/doozit) +1) = V> 

ogr(fix(shoodi/doozit) +1) = gammar; 

oh (fix (shoodi/doozit)+1) = his 


8] 


Rmax, Thbar, 


V, Thbar, 


Thmbar, 


Thmbar i}; 


a a a 








oEs (fix (shoodi/doozit) +1) = Es; 

oe (fix (shoodi/doozit) +1) = e; 

oa (fix (shoodi/doozit) tl) = a; 

ora (fix (shoodi/doozit) +1) = ra; 

erp ( fix(shoodi/dooz2t) +1) = Ep; 

orm tix (snoodi/dooZit) +1) = Tp; 
Thrst(fix(shoodi/doozit)+1l) = Thbar; 
oang (fix (shoodi/doozit) +1) = xbar(3); 


* Increment Number of Orbits 
norbs (fix (shoodi/doozit) +1)=(xbar(3)-x(3))/(2*pi)i 
end 
if rem(shoodi,10*doozit) == 0 
norbs (fix (shoodi/doozit) +1) 


end 


end 
end 


B. VELOCITY SUBROUTINE LISTING 


* NONDIMENTIONALIZED VELOCITY COMPUTATION 

* Establish as Function 

function [Vbarest] = leovel (xbar) 

* Compute Nondimentionalized Velocity 

Vbarest = ((xbar(2))*2+(xbar(1)*xbar(4))%2)%0.9; 


return 


end 


C. DRAG SUBROUTINE LISTING 
~ NONDIMENTIONALIZED DRAG COMPUTATION 


+ Establish as Function 


function [rhogess, Dbarest] = leodrag(xbar, ao, rhoalt, Re, beta, 
Boar, dfact) 


> Determine Density Exponential Term 


rhogess = exp(-((xbar(1) *ao-Re)-rhoalt) /beta) ; 
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Vbar, 


* Compute Nondimentionalized Drag 
Dbarest = (dfact*rhogess*Vbar%2)/(2.0*Bbar) ; 
return 


end 


D. LOCAL CONTRIBUTION SUBROUTINE LISTING 
*£ LOCAL CONTRIBUTION SUBROUTINE 
* Establish as Function 


function [hly -L2,. 13, i4, bo, be; Ly) = leolocal(xbar, mbar, Vbar, 
Thbar, Dbar, alphar) 


* Define Local Computational Variables 


lik = xbar(1)*xbar(4)*xbar(4); 

L2 = 4,0" pi *pi/ (xber (1) *sbar(1) 14 

13 = (Dbar* (xbar(2)/Vbar))/mbar; 

L4 = (Thbar*sin(alphar))/mbar; 

L5 = (2.0*xbar(4)*xbar(2))/xbar(1); 

L6 = (Dbar/(xbar(1)*mbar))*((xbar(1)*xbar(4))/Vbar); 
Li = (Thbar*cos(alphar))/(mbar*xbar(1)); 

return 

end 


E. DIFFERENTIAL VALUES SUBROUTINE LISTING 


2 


POSITIONAL DELTA VALUES ROUTINE 


oe 


Establish as Function 


function [dxl, dx2, dx3, dx4] = leoxdot(xbar, L) 


* Compute Positional Delta Values 


axl = xbar(2); 

ax2 = 1(1) -L(2)-L(3)+L(4) 
dax3 = xbar(4); 

dx4 = -L(5)-L(6)+L(7); 
return 

end 


83 








F. RUNGE-KUTTA SUBROUTINE LISTINGS 
= RUNGE-KUTTA 4TH ORDER COMPUTATION 


* Establish as Function 


function [xbtemp, xbdtemp, nxbar] = leork4a(count, txbar, txbdot, tmane) 


= Perform first part of RK4 


xbtemp(count) = txbar(count); 
xbdtemp (count) = txbdot(count); 
nxbar (count) = xbtemp (count) +0.5*tminc*txbdot (count) ; 


+ RUNGE-KUTTA 4TH ORDER COMPUTATION 
* Establish as Function 
function [xbdtemp, nxbar] = leork4b(count, txbar, txbdot,. xbdot,, tminc) 


+ Perform second part of RK4 


txbdot (count) +2.0*xbdot (count) ; 
txbar (count) +0.5*tminc*xbdot (count) ; 


xbdtemp (count) 
nxbar (count) 


= RUNGE-KUTTA 4TH ORDER COMPUTATION 

* Establish as Function 

function [xbdtemp, nxbar] = leork4c(count, txbar, txbdot, xbdot, tminc) 
® Perform third part of RK4 


txbdot (count)+2.0*xbdot (count); 
txbar (count) +tminc*xbdot (count) ; 


xbdtemp (count) 
nebar (counTc) 


i 


* RUNGE-KUTTA 4TH ORDER COMPUTATION 
= Establish as Function 


function [nxbar] = leork4d(count, txbar, txbdot, xbdot, tminc) 


* Perform fourth part of RK4 


nxbar (count) = txbar (count) ttminc* (xbdot (count) +txbdot (count) )/6; 


G. ORBITAL PARAMETER SUBROUTINE LISTINGS 
+ ORBITAL PARAMETERS 


> Define as Function 
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function[Pr, PV, Pgammar, Ph, PEs] = leoparml(xbar, Tpo, ao, mu) 


* Compute Orbital Parameters 


~ Radius 

Pr = xbar(1)*ao; 

~ Velocity 

PV = (ao* (((xbar(2))*2+(xbar(1)*xbar(4))*2)*0.5))/Tpo; 


* Flight Angle 

Pgammar = atan(xbar(2)/(xbar(1)*xbar(4))); 
* Angular Momentum 

Ph = Pr*PV*cos(Pgammar) ; 

* Specific Energy 

PEs = (((PV%2)/2)-(mu/Pr)); 

return 


end 


* ORBITAL PARAMETERS 

* Define as Function 
function[Pe, Pa, Pra, Prp, PTp, PVimin, PVimax] = leoparm2(Es, h, mu) 
* Compute Orbital Parameters 

* Eccentricity 

etest = 1.0+2.0*Es* ((h/mu)%2); 
Pe = sqrt(etest}); 

~ Semi-major Axis 

Pa = (-mu/(2.0*Es)); 

* Apogee Radius 

Pra = Pa* (irre); 

* Perigee Radius 


Prp = Pa*(1-Pe); 
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+ Period 


PTp = 2.0*pi* (sqrt ((Pa%3) /mu)); 

* Minmum Velocity 

PVimin = ((mu*(1+Pe))/(Pa*(1-Pe)))%0.95; 
+ Maximum Velocity 

PVimax = (mur (l=Pe))/ (Pat (1+Pe)))°0.5; 
return 


end 


H. THRUSTER CONTROL SUBROUTINE LISTINGS 


* THRUSTER FIRING DETERMINATION ROUTINE 
« Establish as Function 
function [Thrust, tfl] = leotflph(Rp, Rmin, 
% Thruster Firing Law 
* Initialize Thrust Logic 
tl = 0; 
*# Turn Thrusters Off 

if Thcheck == Thmax 
& Check Apogee Radius 

if Ra >= Rmax 


Thcheck = 0.0; 
tl = 2; 


end 


end 


~ Turn Thrusters On 
if tl == 0 
if Rp <= Rmin 


Thcheck = Thmax; 
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Rmax, 
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Thmax ) 








end 
end 
Thrust = Thcheck; 
cL. = 217 
return 


end 


= THRUSTER FIRING DETERMINATION ROUTINE 
> Establish as Function 


function [Thrust, tfl] = leotflapn(gammar, grp, V, Thcheck, Thmax, tl, 
Gval) 


* Thruster Firing Law 
~« Turn Thrusters Off 
if tl == 
* Check Flight Angle 
= Turn Thrusters Off 
if gammar >= grp 
Thcheck 


tl = 4;*0 
end 


O02 


end 


*« Turn Thrusters On 


if gammar >=0 
if gammar <= Gval 
Thcheck = Thmax; 
end 
end 


end 
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Thrust = Thcheck; 


trl =-tis 


return 


end 
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